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ABSTRACT 

The speed of sound greatly exceeds typical flow velocities in many stellar and planetary interiors. 
To follow the slow evolution of subsonic motions, various sound-proof models attempt to remove 
fast acoustic waves whilst retaining stratified convection and buoyancy dynamics. In astrophysics, 
anelastic models typically receive the most attention in the class of sound-filtered stratified models. 
Generally, anelastic models remain valid in nearly adiabatically stratified regions like stellar convec- 
tion zones, but may break down in strongly sub-adiabatic, stably stratified layers common in stellar 
radiative zones. However, studying stellar rotation, circulation, and dynamos requires understanding 
the complex coupling between convection and radiative zones, and this re quires robust equatio ns valid 
in both regimes. Here we extend the analysis of equation sets begun in | Brown et all (|2012[ ). which 
studied anelastic models, to two types of pseudo-incompressible models. This class of models has 
received attention in atmospheric applications, and more recently in studies of white-dwarf super- 
novae progenitors. We demonstrate that one model conserves energy but the other does not. We use 
Lagrangian variational methods to extend the energy conserving model to a general equation of state, 
and dub the resulting equation set the Generalized Pseudo-Incompressible (GPI) model. We show 
that the GPI equations suitably capture low frequency phenomena in both convection and radiative 
zones in stars and other stratified systems, and we provide recommendations for converting low-Mach 
number codes to this equation set. 
Keywords: stars:interiors - Sununterior 



1. INTRODUCTION AND MOTIVATION 

In astrophysical fluid dynamics, the relevant timcscalcs 
are often substantially longer than the sound crossing 
time of the system. This particularly holds true for con- 
vection deep in stellar interiors where the flows are very 
subsonic. Near the base of the solar convection zone the 
sound speed is about 220 km/s, while the convective ve- 
locities are likely of order hundreds of meters per second. 
Following the evolution of sound directly imposes crip- 
pling computational limits on simulations of such flows, 
as their evolution times are typically many convective 
turnover times, each of which is often several thousand 
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sound times. 

For numerical stability, an explicit time-integration 
scheme must satisfy the Courant-Friedrichs-Lewy (CFL) 
condition, that the time-step must be smaller than 
the shortest timescale admitted by the equations. For 
the non-rotating Euler equations in a stably stratified 
medium, this is 



At < 



Ax Aa; 1 

c ' u ' N 



(1) 



where At represents the time-step size, Aa; is the small- 
est resolved length scale in the calculation, and cq, uq, 
and Nq represent characteristic sound speed, flow veloc- 
ity, and Brunt- Valsala frequency respectively. The Mach 
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number gives the ratio of the first two timescales in equa- 
tion dTJ 



Ma 



Co 



(2) 



which is typically very small in stellar interiors. If we 
make the order of magnitude estimates c\ ~ gH and 
Nq ~ g/H, where g, and H represent gravitational ac- 
celeration and density scale height respectively, then the 
ratio of the first and third timescales in equation ([1]) is 



Aa; 



(3) 



which is also small in high-resolution simulations of stel- 
lar interiors. 

So called "sound-proof" models address this separa- 
tion of time scales by starting with the Euler equations 
and filtering out fast, high-frequency sound waves while 
retaining compressible motions on slower time scales due 
to gravitational stratification. These motions include 
gravity waves in stably stratified regions and convection 
in unstably stratified regions. The CFL condition for 
sound-proof equations only requires 



. Ax 1 

At <- mm < , — 

u Nq 



(4) 



which much-less severely restricts efficiency than equa- 
tion ([1]) for the fully compressible Euler equations. In as- 
trophysical and geophysical settings, the most commonly 
emplo yed "sound-proof " models are the anelast i c equa- 
tions (jBatchelori fl953t IQgura fc PhilhfJsl H96l IGoughl 
1 1969T ) . Fundamentally, anelastic models filter sound 
waves by modifying the continuity equation of the fully 
compressible Euler equations. Formally, this anelastic 
approximation is onl y valid for an a diabatic or nearly adi- 
abatk atmosphere (jGoughl 119691 : iBraginskv fc Roberts! 
[19951) . 

The pseudo-incompressible models — with a modified 
pressure, rather than density, equation — give an al- 
ternate app roach to sound proofing the Navier-Stokes 
equations. IDurranI (|1989[) first proposed this class of 
models, which the astrophysical fluid dynamics com- 
munity mo re recently adapted to a number of applica - 
tions (e.g., lAlmgren et afll2006albl iZingale et afll2009l ). 
A pseudo-incompressi ble model finds par ticular use in 
the MAESTRO code (jNonaka et al.ll2010h . The atmo- 
spheric sciences community has also extensively explored 
the properties of gravity waves and stable-layer dynamics 
in pseudo-incompressible models, with several compar- 
isons against the properties of anelastic models (IDurranI 
1989L 12008^ INance fc Durrani [Tool lAchatz et al.l 120101 : 
Klein et al 



I2010h 



In iBrown et al.l (|2012f) (hereinafter Part I) we dis- 
cussed the energy conserving properties of various anelas- 
tic models in stably stratified, sub-adiabatic atmospheres 
such as are found in stellar radiative zones. We found 
that some widely used anelastic models violate energy 
conservation. This behaviour escaped attention previ- 
ously as in bounded atmospheres, the norm for most sim- 
ulations, energy-violating anelastic models instead con- 
serve a pseudo-energy (an energy-like quadratic invariant 
with incorrect stratification weighting) . Internal gravity 



wave eigcnfunctions in energy- violating anelastic mod- 
els can differ significantly from the fully compressible 
results, and neither energy nor pseudo-energy remain 
conserved when nonlinear dynamics become important. 
Anelastic models that correctly conserve energy have 
modified momentum equations, a theme which reoccurs 
here in our study of pseudo-incompressible models. 

Here we explore two pseudo-incompressible models 
(Sj2|), one used in the atmospheric community (PI equa- 
tions) , and a new variant being used in the astrophysical 
community that includes more general equations of state 
(LM equations). We find that the PI equations conserve 
energy but that the LM equations do not. We addi- 
tionally derive a new subsonic flow model based on a 
constrained Lagrangian analysis of the full compressible 
flow (Sj3] & SJ). We call these the Generalized Pseudo- 
Incompressible (GPI) equations. The GPI equations 
conserves energy and correctly generalizes the pseudo- 
incompressible approximation for an arbitrary equation 
of state including astrophysically relevant situations such 
as radiation hydrodynamics. After analyzing the general 
properties of the GPI, PI, and LM equations, we special- 
ize to the case of an ideal gas equation of state to explore 
the behavior of these differing models in bounded at- 
mospheres (JJ5J) and perform numerical simulations that 
show the difference between the LM and GPI equations 
(Sj6]). The implications of these findings for simulations of 
stellar interiors is discussed in Sj7] The essential results of 
this paper are the derivation of the Eulcr-Lagrange equa- 
tion (|73|) for any general non-dissipative fluid flow, and 
application of equation (|73[) to derive the Generalized 
Psuedo-Incompressible equations as expressed in equa- 
tions (|99l) - (|101[) . The reader who is primarily interested 
in implementing energy-conserving subsonic models with 
a general equation of state should read [H [6] and [7] 

2. BACKGROUND AND MODEL EQUATIONS 

2.1. Thermodynamics and Stratification 

We begin our discussion with thermodynamics and the 
geometric conservation laws of mass and entropy 

dp 

dt 
ds 

di 



V • (pu) = 



u • V,s 



0. 



(5) 
(6) 



where u(t,x) denotes the Eulerian (fixed spatial coor- 
dinates) fluid velocity, p(t, x) denotes the density, and 
s(t,x) represents the specific entropy. The main result 
of this paper provides a framework for producing appro- 
priate evolution equations for the fluid velocity under 
different physical assumptions and approximations. At 
this stage, u(t, x) represents an arbitrary flow field. 

The combined First and Second Laws of Thermody- 
namics relate small changes in density, entropy and spe- 
cific internal energy, e, via, 



Tds 



■dp, 



(7) 



Equation ([7]) suggests that the internal energy naturally 
depends on the density and entropy. An equation of 
state closes the thermodynamic description of the sys- 
tem, which comes in the form of an algebraic (instan- 
taneously valid at every point in space) relation be- 
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tween any three thermodynamic variables; for example 
p = p(p,s), or more fundamentally e = e(p,s). There- 
fore, 



(8) 



define the pressure, p, and temperature, T. We consider 
the two partial derivatives of internal energy e at con- 
stant entropy s and density p respectively. 
We also find it useful to define the specific enthalpy, 



h 



P 
P 



dh 



Tds + ^. 
P 



(9) 



The differential form of equation © suggests that the 
internal enthalpy depends naturally on the pressure and 
entropy. 

If using h(p, s), rather than e(p, s), and assuming {p, s} 
evolve according to equations (JSJ & ([6]), then determining 
the pressure requires implicitly solving 



1 _ f dh(p,s) 
P V dp 



(10) 



While the distinction between internal energy and en- 
thalpy appears like a mere reorganization of the ther- 
modynamic variables, using the enthalpy and pressure 
significantly aids in interpreting low Mach number ap- 
proximations to dynamical models. 

Finally, in addition to equations (jSJ & (|10|) . equa- 
tions ©-(J?]) also imply an explicit local pressure evo- 
lution equation. In differential form 



dp 



IK 



(11) 



Therefore, interpreting the differentials as convective 
time derivatives produces, 

d_ 

Of 



TT+u-V )p + T 1 (p,p)pV -u = 0, (12) 



where 



p \dp 



(13) 



specifies the first adiabatic exponent. The adiabatic 
sound speed also follows from equation (fTS]) . 



c 2 {p, p) = 



dp 
dp~ 



= Ti(P,p) 



P 



(14) 



In Part I, we restricted our analysis to a monatomic 
ideal gas with Ti = 7 = 5/3. Here, we derive a general 
sound-proof hydrodynamic model with a general equa- 
tion of state. Using general equations of state provides 
a two-fold advantage over the ideal gas model. First, 
many astrophysical application simply require more com- 
plex physics than ideal gases can model. In Appendix [X] 
we give an example of an astrophysically relevant equa- 
tion of state incorporating an ideal gas and blackbody 
radiation, with the general form for Ti given in equa- 
tion (|A4[) . This is of particular interest when studying 
convection in massive stars, where near-Eddington lumi- 
nosities lead to strongly interacting mixtures of radiation 



and matter, but the fluid dynam ics remain very subsonic 
(|Cantiello Ik Braithwaitel l201li ). One can find these and 
other examples, including partially ionized gasses criti- 
cal fo r stellar interiors, i n textbooks on stellar structure 
(e.g.. ICox fc Giulil 119691 . Chapter 9). Second, keeping 
a general equation of state helps clarify and unify the 
thermodynamic interpretation of a number of concepts 
that arise when considering the energetics of waves and 
instabilities in stratified atmospheres. We also show in 
Appendix [Al how to implement the GPI equations ([99)) 
(I101D using this more complex equation of state. 

One of the main themes of this paper and Part I fo- 
cuses on the consequences of gravitational stratification 
for low Mach number dynamics. Therefore, we intro- 
duce a hydrostatically balanced, stratified reference at- 
mosphere with background density pa , pressure po , tem- 
perature To, and entropy sq that vary only in the direc- 
tion of gravity 



Vp = -pa V(; 



(15) 



The potential function <p(x) determines the local gravi- 
tational acceleration g = — V0. 

Part I discusses anelastic models, which eliminate 
acoustic modes from the dynamics by employing the den- 
sity constraint 



\p~Po\ < Po 



V • (poti) = 0. 



(16) 



Equation (|16|) rests on the premise that dynamical den- 
sity fluctuations remain small compared to the back- 
ground. This assumption indeed filters sound modes, but 
is only formally justified for nearly adiabatically strati- 
fied atmospheres. 

Beginning with iDurranl (|1989f ) , a number of authors 
have recognized that the pressure field more naturally 
distinguishes between acoustic and low Mach number dy- 
namics. The pressure field provides the restoring force 
for acoustic oscillations, and this implies that small pres- 
sure fluctuations, 



(17) 



characterize low Mach number dynamics on small scales, 
rather than small density fluctuations as in equation (|16| . 

Therefore, replacing the pressure field with the static 
background reference field, the pseudo-incompressible 
models in this paper replace the pressure evolution equa- 
tion (I12p with a velocity divergence equation 



V • u 



u • Vp 

Ti(po,p)po 



0. 



(18) 



Moreover, for the case of constant T± — 7, we may 
rewrite equation (|18|) into a form similar to equation (|1( 



where 



v • (A>«) = 0, 



(19) 



(20) 



We consider equation (|18| the defining characteristic of 
pseudo-incompressible models — as opposed to anelastic 
models. In adiabatically stratified atmospheres, sq — 
constant implies po oc pg , and 



V • (Ajti) oc V • (p u) = 0. 



(21) 
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In these atmospheres, the pseudo-incompressible con- 
straint equation (|18l) reduces to the more common anelas- 
tic constraint, which we studied in Part I. 

Here we consider two different pseudo-incompressible 
models. The different notation and different thermody- 
namics used in the various treatments can lead to ambi- 
guity over the equivalence or differences between models. 
To avoid this, as in Part I, we write each model using as 
consistent a notation as possible. Practical numerical 
or computational differences can arise when solving dif- 
ferent transformations of the same fundamental model, 
but these issues mostly lie beyond our current scope. 
Therefore, we consider two models identical if one can 
bring them into the same form by legitimate mathemat- 
ical transformation (without approximation), and/or by 
a possible change of notation. 

2.2. Fully compressible Euler equations 

The thermodynamic equations ([5"])- (|12[) apply with re- 
spect to an arbitrary flow field, u(t,x). For fully com- 
pressible dynamics neglecting dissipation and sources of 
heat, the fully compressible Euler (FC) equations govern 
the flow evolution. 



( ^ + u ■ V 



= 



dp 
dt 
ds 
di 



V • (pu) = 



u • Vs 



0. 



(22) 

(23) 
(24) 



An equation of state p = p(p, s) closes equations (1221 
(|24p . In this paper, with conservation of energy a pri- 
mary theme, we limit our analysis to strictly isentropic 
flow, ds = (equation (|24jl ). This allows us to ensure 
that no mechanical violations of energy balance occur. 
Our more general interest lies in proper energy budgets 
for non-iscntropic flows, but including dissipative terms 
would complicate the analysis presented here, and we 
intend to consider such effects in future work. 

Equations (f2"2"T) -(f24" |) imply the conservation of total en- 
ergy 



dE 

— + V-[u(E + p)} = 



where 



E 



p\u\ 



+ P e (p, s ) + P<t>- 



(25) 



(26) 



The FC equations conserve energy with an arbitrary 
equation of state, assuming dt<p = 0. 

The form of the potential energy in equa- 
tions ([2"5]) & (Unj) leads to some confusion when 
comparing linear and nonlinear dynamics. For nonlinear 
states, the potential energy appears to determine very 
little about how the system may or may not extract 
energy from a background stratification and convert it 
into fluid motion. Alternatively, by perturbing about a 
background hydrostatic stratification, the equations of 
motion possess a well-known energy-conservation prin- 
ciple that adopts a quadratic form in the perturbation 
variables. 



We reconcile these facts by noting that equa- 
tions J5]) & (|6]) together imply 



dt 



V • [puf(s)} = 0, 



(27) 



for any arbitrary function exclusively depending on en- 
tropy, f(s). 

Equation (|2"T|) implies the existence of a family of con- 
served free energies, 



F = E - pf(s) 



p\u\ 



+ p[e(p,s)- f(s) + </>], (28) 



where we subtract an arbitrary "ground-state" energy, 
represented by f(s). 

Given a hydrostatically balanced background stratifi- 
cation {p (x),po(x), so(x)}, we may choose f(s) such 
that the linear contribution from perturbations to the 
free energy cancel identically. The character of the 
quadratic terms then determines the linear stability, and 
possibly even the nonlinear stability. 

The i dea of available potential energy goes back to 
Lorenzl ( 19551) in the case of a Boussinesq fluid, and 
Andrews! (|1981[ ) computed the quantity for a compress- 
ible system by analyzing the dynamical equations di- 
rectly. Using the potential energy, we show in Ap- 
pendix [B] that the functional relations 

f(s (x)) - h (x) + 4>{x), f'(s (x)) = T (x), (29) 

define the appropriate choice of ground state energy den- 
sity required to cancel contributions from linear pertur- 
bations. Mathematically, ho + (j>, represents the "pull- 
back" of / with respect to So- We may use equation (|29[) 
to find a global function, /(s), if sq is a one-to-one func- 
tion of the gravitational potential, (f>. 

Equation ([29]) renders F quadratic to lowest-order in 
the perturbations around background state. That is, 



p\u\ 



+ A-po(x), 



(30) 



where A represents the available potential energy. To 
leading order, 



A = p (e(p, s) + 4>{x) - f{s)) + p (x) 



{P-Po? t poN§(s-s ) 2 
2/OqCo 



where, 



N* 



2 Vs 



Po 



h.o.t., 



(31) 
(32) 

(33) 



defines the background Brunt- Vaisala frequency, and cq 
represents the background sound speed. For linear dy- 
namics, we may neglect the higher-order terms in equa- 
tion pip , and obtain the well-known quadratic energy 
invariant in equation (|32[) . In this case, we see clearly 
that iVg > implies F > — po, and the system cannot 
create kinetic energy indefinitely. We also see that, to 
leading order, pressure perturbations always produce a 
positive contribution to the free energy. 

Additionally, we show in Appendix [C] that one may 
rewrite equation (|33l) into the alternative general form 



(34) 
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where c p gives the specific heat capacity at constant pres- 
sure, and ut gives the thermal expansion coefficient at 
constant pressure, 



t(— 
\bt 



CUT 



dp 
df 



(35) 



Equation (jM]) relates entropy gradients to the potential 
density gradient in equation (|33p . and simplifies to the 
more typical expression in the case of an ideal gas where, 
Tax = 1, and c p is constant. 

Along with energy, entropy and mass, non-dissipative 
flow dynamics conserves the Ertel potential vorticity 



such that, 



at 



(36) 



(37) 



Equation (|37[) represents a fundamentally important re- 
sult in continuum fluid mechanics and constrains many 
important processes. Conservation of potential vortic- 
ity finds many useful applicati ons across geophy s ical and 
astrophysical fluid dynamics ( jBrethertonl I197G : iSalmonl 
Il988f ) and neglecting equation (|3"T|) admits similar errors 
as neglecting energy conservation in equation (I25|) . 

2.3. Pseudo-incompressible equations 

iDurranl (jl989t ) first proposed pseudo-incompressible 
models as a response to inadequacies of anelastic models. 
As Part I discusses in detail, anelastic models faithfully 
describe nearly adiabatic convection zones whilst filtering 
fast acoustic motions from the dynamics. Additionally, 
different anelastic models perform with differing results 
in stably stratified zones. 

The anelastic equations filter sound waves with a con- 
dition equivalent to requiring the density field to in- 
stantaneously match the background stratification. As 
we discuss in the introduction, however, local pressure 
equilibrium proves the more physically well-motivated 
assumption when wishing to sound-proof a dynamical 
model. The distinction between density and pressure 
in approximations becomes especially clear for atmo- 
spheres with significant sub-adiabatic (stable to convec- 
tion) stratification. In the Sun, this particular diffi- 
culty arises in modeling the transition from the con- 
vective outer envelope to the stable radiative interior. 
This region of penetration and overshoot, known as the 
tachoclinc, likely plays an important role in the so- 
lar dynamo. Understanding the coupled dynamics be- 
tween convection and gravity waves in this region is vi- 
tally important. Sufficiently deeply in the solar radia- 
tive zone, the assumption of near adiabaticity underly- 
ing the anelastic approximation eventually breaks down, 
even though the flow remains in a significantly low Mach 
number regime. 

The work of IDurranl (fl989h sought to overcome the 
difficulty of transitions to strong stable stratifications 
when modeling similar situations in earth's troposphere 
and stratosphere. The main innovation in that work 
assumed rapid pressure, rather than density, equilibra- 
tion. However, even though this model, dubbed the 



pseudo-incompressible (PI) equations, significantly gen- 
eralize the anelastic equations, it still implicitly assumes 
an ideal gas equation of state. This assumption leaves 
open the question of how to generalize the idea of pres- 
sure balance to the more complex situations encountered 
in many astrophysical applications. In this section, we 
outline the basic assumptions of the PI equations with 
an ideal equation of state. After establishing the ba- 
sic idea, we use a general variational framework (Q to 
move beyond simple equations of state whilst retaining 
the appealing features of the PI equations in 

The FC equations from ^2.2l rcauire adjustment to cor- 
rectly accommodate the pressure-restricted divergence 
equation (|19|) . The PI equations evolves according to 



+ (p-po)9, 



g + V-(p«) =0 

v-(Aju) = o 



(38) 

(39) 
(40) 

where D/Dt = dt + u • V, and p' = p — pa repre- 
sents the pressure fluctuations around a given hydro- 
static background. We recall the definition (3q = p^ 1 
(equation (|2"0"1) ). Equations (f3"B")) - (|4"0")) assume an ideal gas 
with constant ratio of specific heats 7. Ostensibly, the 
pressure-gradient term in the momentum equation (|38|) 
takes an unusual form when c o mpar ed to the FC mo- 
mentum equation ([22]) . iDur ran (1989) originally derived 
equation (|38p by writing the FC momentum equation in 
terms of the Exner function and potential temperature 
respectively, 



P 

Prof 



T 

7T 



(41) 



where p rc f is a constant reference pressure. The only 
approximation comes from assuming small Exner func- 
tion fluctuations in the pressure evolution equation. Q 
As we show in Jj4j the treatment of pressure in momen- 
tum equation (|38|) provides exactly the correct behavior 
given the pseudo- incompressible divergence constraint. 
Equations (|38|) - (|40|) also conserve entropy, assuming an 
equation of state of the form s = s(po, p). The most im- 
portant physical assumption of the PI equations is that 
they ignore pressure fluctuations everywhere except the 
momentum equation. 

The conservation of energy gives one important self- 
consistent feature of the PI equations, 



BE 



pi 



at 



v 



where, 



u [Epi + p' 



p\u\ 



7 



7-1 



Po 



Epi = 



= 0, (42) 



(43) 



We note that 7Po/(7 — 1) = Po n o coincides with the back- 
ground enthalpy density for an ideal gas. Equation (j42|) 
states that filtering acoustic motions from the dynam- 
ics should not produce energy balance violations in low 



1 Tran s form ing the equations in IDurranl l)1989D into equa- 
tions j2H)-gQj 

requires identifying p* with the mass density p. 
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Mach number flow. Wc expand on the physical reasons 
for this important fact in 21 

Similar to the FC equations, the PI equations also con- 
serve a free energy that assumes a quadratic form for 
small perturbations, 



F = pa- 



u 



u 



■Po- 



Po n q{s ~ sp)" 2 

2|Vs | 2 
9 2 (P-Po) 2 



ipo 



h.o.t. 



h.o.t. 



(44) 



where we only retain the leading-order quadratic form for 
linear wave dynamics. Importantly, equation (|44[) only 
omits free potential energy associated with acoustic oscil- 
lations. Dropping pressure dependencies allows us to use 
only the density perturbations, for sources of potential 
energy. Equation ()29[) determines the adjustment, /(s), 
to equation (|4"3"]l that produces equation (|4"4")) . just as for 
the FC equations. The conservation of equation (|44[) for 
small-amplitude gravity waves makes an important point 
of comparison to other models in iJBJ 

In addition to properly conserving energy, the PI equa- 
tions also consistently account for potential vorticity in 
the reduced dynamics. In this situation, the entropy gra- 
dient reflects the fact that p rts po, 



Vpo 

7Po 



Vp 

P 



and 



d 

di + U 



Vs'Vxii 

P 



= 0. 



(45) 



(46) 



2.4. MAESTRO low Mach number equations 

In the last few years, a new subsonic model came 
into the astrophysical community from work modeling 
Type-la supcrnovac progenitors. These low-Mach num- 
ber (LM) equations , resemble the PI equations but 
with the important feature of a generalized equation 
of state needed to model the complex thermodynam- 
ics of converting whi te- dwarf interiors (j Almgrenl 120001 : 
lAlmgren et all [2006al lbl. l2008h . The MAESTRO code 
implements the LM equations in an adaptive resolution 
framewor k suited to the needs of multiscale convection 
problems (iZingale et al.|[2009l: iNmlaka et al .1 120 lOl [20121: 
iMalone et al.ll2011[) . 

The LM equations use the same momentum and 
density equations as the standard-form Euler equa- 
tions ((SJ) & (|22|) . However, the velocity-divergence equa- 
tion (fT5|) replaces the pressure equation (TT21 . Wc write 
the LM equations set as 



P-q£ = -Yp + (p- pa)g 

dp 
dt 



V • u 



(pu) = 
Li(po,p)po 



= 0, 



(47) 
(48) 
(49) 



where we again use p' to denote the pressure perturba- 
tion. We here exclude many dissipativc and multi-species 
effects that are important in the modeling of dynamics 
in actual white-dwarf interiors; including these effects 



would not change the fundamental conclusions we draw 
here about the LM equations. 

The LM momentum equation (|47[) is motivated by scal- 
ing the FC momentum equation (written in terms of p 
and p) with respect to the characteristic velocity U le t, 
the timescale i re f = U re i/g, and th e length scale L re f = 
J7 ref /t ref = U? e{ /g (|Almgren et alJl2006aft . The density 
and pressure scale with respect to p re f and p re f , and the 
characteristic pressure scale height H rc { = p re f/(Pref <?) 
is defined such that L TC f/H rc { = 0Ma 2 . These choices 
imply that the non-dimensional momentum equation be- 
comes 



Du* 
~DtI 



1 



Ma- 



■V*(p* -po,*) + (p* - P*,o)0», (50) 



where the * subscripts represent non-dimensional scaled 
quantities. The non-dimensional pressure fluctuations 
must remain 0(Ma 2 ) in order to balance time evolution 
and buoyancy, and omitting the evolution of these fluc- 
tuations from equation (|12|) filters acoustic modes and 
allows for volume cha nges by moving throu gh the back- 
ground stratification (|Almgren et al.ll2006al) . However, 
equation ( 50 ) gives a different momentum equation than 
equation (|38[) because transforming the original thermo- 



dynamic variables in the FC equations, and then approx- 
imating the result, gives a different result than making 
the approximations first and then transforming to new 
variables. This particular derivation neglects possible 
changes in buoyancy due to momentum exchange — but 
not energy exchange — between the low Mach number 
flow and rapid sound waves. 

In trying to derive an energy conservation principle for 
the LM equations, the following form gives the closest 
possible analogue to equations ([2"B"]) & (|4*2"|) . 



dE 



LM 



dt 



V ■ (m£ L m) = -u-Vp', 



(51) 



where 



E 



LM 



p\U\ 



+ pcf> + ph(p , s) = E P1 + ph(p„, s), (52) 



and h(p , s) evolves according to equation ©. The term 
on the right-hand side of equation (|5ip does not reduce 
to an exact divergence, and represents an uncontrolled 
energy source or sink. For an ideal gas with constant 
ratio of specific heats, 7, equation (|5T1) becomes 



dE 



pi 



dt 



+ v 



u E 



'j-pi 



7 



1 



Po 



= — u • Vp', (53) 



which only differs from equation (|42[) by p' V • u on the 
right-hand side. For general flows, the LM equations do 
not conserve any energy E. 

Along with energy, the LM equations also do not con- 
serve potential vorticity. To show this, we use equa- 
tion (|C5j) 

ds dp dp 

la T — = , 

c p Tip p 

where ut is the coefficient of thermal expansion given in 
equation (|55|) . Because of the pressure constraint, p k, 
Po, the entropy for the LM equations only depends on 
density p and po, thus s = s(po(x), p). The entropy s 
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still satisfies equation ©. Therefore, q = Vs-(V Xu)/p 
still defines the potential vorticity, using 



Vs 



Tcht \ r iPo P / 



(54) 



where c p ,T, or, and Ti all implicitly depend on po and 
p. Therefore, 

9g 



9i ' ^ p 3 raTriPo 



Vpo • (Vp X Vp') . (55) 



If we assume an ideal gas equation of state, then equa- 
tion ([54"]l reduces to equation (|45|) and equation (|55l) re- 
duces to 



^9 



P 3 Po 



Vpo • (Vp X Vp') . (56) 



No way exists to eliminate the right-hand side of equa- 
tions ([55]) or ([56| . The LM equations do not conserve 
potential vorticity for general motions. 

Problems remain even in the linear regime. The right- 
hand side of equation ([55]) docs vanish if we consider only 
linear perturbations, so linear perturbations do conserve 
potential vorticity, which vanishes for gravity waves. 
They do not however conserve energy. In the linear 
regime, the conservation of s and the definition of the 
buoyancy frequency in equation (|34p implies 



dE 
~9t 



V • (up ) = — p - 



u ■ Vp 



(57) 



p a Ti(p ,p a )' 

where E follows from equation (|44|) just as in the PI 
equations. 

(58) 
(59) 

(60) 



E = K 
Po\u 



2 



K 



U = 



9 2 (P-PQ? 
2N%p Q 



where we use the density, rather than entropy, form of 
equation (|44|) . Even for linear motions the right-hand 
side of equation (|57[) remains non-zero. 

Although the LM equations do not conserve energy, 
the linearized equations do possess a quadratic invari- 
ant, that we call the pseudo-energy. As in Part I, the 
presence of a pseudo-energy can produce misleading ef- 
fects since it allows bounded solutions in some dynamical 
regimes. However, bounded solutions do not necessarily 
imply accurate solutions, and as we found in Part I, these 
pseudo-energy solutions may contain significant errors in 
many situations. We now derive the conserved pseudo- 
energy. 

Assuming po and po only depend on the directio n of 
g, which we call z, we follow lAlmgren e~a l. (2006a|) and 
define a more general 



log A) 



dpo 



i 9 -f- m 



p Ti(p ,p ) j u 

For the simple case of an ideal gas with constant 7, equa- 
tion (fBTj) recovers f3o = p^ 1 . Using the general form of 
/3o in equation (|57|) produces 

dPE 



dt 



V • (up'Po) = 0, 



(62) 



where PE represents the quadratic invariant pseudo- 
energy 



PE = (3 Q E = PK + PU, 
PK = p Q K= PoP °W\ 



PU = (3 U = 



9 2 Mp-po) 2 
2A 2 Po 



(63) 
(64) 

(65) 



Equation (f6"2")l implies that linear motions in the LM 
equations conserve pseudo-energy PE. 

The pseudo-energy and energy only coincide for con- 
stant /So, which implies an atmosphere with no back- 
ground stratification. For an adiabatic atmosphere /3q cx 
po, and the LM do not reduce to anelastic equations. 
In SJHI we numerically calculate the evolution of the en- 
ergies (equations (|58 |) --(|60 |) ) and pseudo-energies (equa- 
tions (f6"3")) l|65p) for a simple test problem with both the 
PI and LM equations. 

3. LAGRANGIAN ANALYSIS 

The analysis in Sj2] makes clear that the form of the 
momentum equation determines if a given sound-proof 
model will correctly conserve energy. Generally, as we 
found in Part I for anelastic models, sound-proof equa- 
tions using an unmodified Euler momentum equation do 
not conserve energy. In this section we develop the tools 
which allow one to correctly and consistently determine 
the momentum equation for a soundproof model with 
general equations of state. 

3.1. Eulerian Action Principle 

Students in fluid mechanics commonly first learn to 
derive equation (f22|) by appea ling to Newton's laws of 
motion for a continuum media (Landau &: Lifshitzlll959l : 
iKundu. Cohen fc Dowlind |2012[ ). This type of analysis 
typically considers fluid elements either moving with the 
flow (Lagrangian formulation), or in fixed volumes (Eule- 
rian formulation) and analyzes the forces acting on these 
fluid elements. Given the analysis of forces, typically 
both of these approaches fall under the scope of Newto- 
nian mechanics, regardless of where one places the spatial 
coordinates. Alternatively, Lagrangian mechanics relies 
on Hamilton's principle of stationary action to reformu- 
late a dynamical equations of motion without recourse 
to specific forces. Here wc derive the fully compressible 
Euler momentum equation (|22[) using Lagrangian me- 
chanics. 

In contrast to Newtonian dynamics, Lagrangian me- 
chanics defines a Lagrangian function £, which is related 
to the energy at each point of the medium. The action, 
<S, is the space-time integral of the Lagrangian density. 
Hamilton's principle states that the equation of motion is 
equivalent to the stationarity of the action with respect 
to the variables of the system. 

Hamilton's principle assumes a simple form and inter- 
pretation when considering particle mechanics. The sit- 
uation in continuum fluid mechanics is somewhat more 
complicated. In our case, the problem corresponds to the 
stationarity of the action 



C(u, p, s) d 3 x dt, 



(66) 
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such that the density and entropy satisfy equa- 
tions ([5|) & (j6)). We leave the particular form of the 
Lagrangian unspecified until the next section. Given the 
similarities with particle mechanics, the Lagrangian co- 
ordinate formulation of a fluid lends its elf well to solution 
by Hamilton's principle (jSalmonlll988f ). In this case, af- 
ter extremizing the action, one must take the trouble of 
recasting the dynamic equations into an fixed coordinate 
(Eulerian) form to compare with the Newtonian-derived 
equation (|22|). Alternatively, one may consider an Eule- 
rian version of Hamilton's action principle on a fixed spa- 
tial coordinate system. In this case, the added difficulty 
results from properly enforcing the geometric conserva- 
tion constraints of mass and entropy equations ([5]) & (|6]). 

In Eulerian coordinates, the volume measure remains 
fixed. Therefore, Hamilton's principle becomes 



SS 



dC 

du 



• Su 



dC 
dp 



Sp ■ 



dC 



Ss d 3 a;dt = 0. (67) 



Two constraints relate the virtual displacements in equa- 
tion ([67]) . 



^ + V • {pSu + uSp) = 

— - + u ■ VSs + Su ■ Vs = 0. 

dt 



(68) 
(69) 



Equations ([68")) & (|69|) follow by linearizing equa- 
tions (O & ([HI) around a given phase-space trajectory 
{u(t, x),p(t, x), s(t, x)}. Two main methods exist for en- 
forcing equations (|68|) & (|69[) in Eulerian variables. The 
first method uses Lagrange multipliers to constrain the 
action directly. This method of Lin constraints, requires 
the added complicat ion of Euler (a.k.a. Clebs ch) poten- 
tials for the velocity (|Cendra fc Mars dcn 1987]). The sec- 
ond method for enforcing equations JS]) & © amounts 
to directly parameterizing equations (|B"8|) & (f69"|) . A well- 
known result from linear analysis, the following pertur- 
bations 



Su = 


d£ 

-£+u-V£-£-Vu 

ot 


(70) 


Sp = 


-v • {pi) 


(71) 


Ss = 


-£ ■ Vs 


(72) 



directly solve equations (f68|) & (|69l) for any displace- 
ment £, and any flow fi eld it assuming p and s sat- 
isfy equations © & © (jEckartl fl960t INewcombl Il96l 
Kulsrud 2005). The three-dimensional infinitesimal dis- 
placement, £, generates the five-dimensional constrained 
phase space. One may easily confirm equations (f70")) -(f72" j) 
if £ = uSt, Su = (d t u)St, Sp = (d t p)St, and Ss = (d t s)St. 
For this special case, equation ((70)) reduces to a triv- 
ial statement, and equations (f7"T|) & (f72j) become equa- 
tions ([5]) & (JU) respectively. 

Using equations (|T0"|) - (f72"j) in equation (|6T|) and inte- 
grating by-parts to isolate £ gives the Euler-Lagrangc 
equations for an ideal fluid in a fixed spatial Cartesian 
coordinate system, 



dfdC\ d f dC 
dt\dui J dxj \ 3 dui 
d fdC\ fdC\ ds 



dxi \ dp 



ds J dxi 



dC \ duj 
duj ) dxt 

0, 



(73) 



where we use Einstein notation and sum over repeated 
subscripts. 

Equation (|73[) represents the Euler-Lagrange equa- 
tions of motion for a general non-dissipative fluid. At 
this stage, we specify nothing about the specific form 
of the Lagrangian density function. This master equa- 
tion allows us to derive the FC momentum equation (|22p 
from an alternative perspective, but more importantly 
allows us to make approximations in a general way and 
distill dynamical models that remain consistent with re- 
gard to energy, and potential vorticity budgeting, as we 
now demonstrate. 



3.2. Energy conservation 
We define the momentum density 

dC 
du 



3 = 



(74) 



With this, conservation of energy follows by contracting 
equation ([73| with the velocity field, u, 



dj _ 

d(u -j-C) 

dt 



[u (u ■ j)} 



pU - V d~p 

dC 

[U - J - p dp- 



_ dL 

ds 



= 0, 



(75) 



where the last line assumes that the Lagrangian C does 
not depend on time explicitly. We recognize the Lcgendrc 
transforms of the Lagrangian density with respect to the 
velocity and density give the "Hamiltonian" Ji, and the 
"generalized pressure" V, respectively, 

dC 

H = u-j-C, V = C-p—. (76) 

dp 

From equations (|75|) & (J75J), the general statement of 
local energy conservation follows in conservative form, 

dU 
dt 

The primary utility of equations (|T3"|) & ([77]) lies in the 
unspecified form of the Lagrangian density C. If we 
choose the traditional form for C (equation ([86]) ), as we 
will in §3.4[ then equation ([77]) leads to equation (|25[) 
and we have energy conservation. Furthermore, just as 
in §2.21 and equation (|27|). we may define the conserved 
free energy 

F = U-pf{s) (78) 
Assuming the following relationship 



+ V ■ [u {% + V)} = 0. 



(77) 



for any f(s) 



/(*(*)) - f 



(79) 



u=0,p=p (x),s=s (x) 



eliminates all linear contributions to the available po- 
tential energy, thus generalizing equation (|2T)1) . We also 
encourage the reader to compare the Lagrangian formu- 
lati on for fluid mech anics with the Hamiltonian formula- 
tion [Morriioi (|1998t) . 

3.3. Potential vorticity conservation 

We now turn to the conservation of potential vorticity. 
We begin by noting that for any tp, if 

£ = ^2iZ£, (80) 
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then Sp — Ss = 0, but Su ^ 0. 

The interpretation of potential vorticity becomes clear 
when one realizes that equations (|5|) & ([6]) imply that 
density and entropy act as dynamic coordinate vari- 
ables. Equations (fT0")) -(fT2" j) show that d t $, generates 
three-dimensional virtual displacement velocity, and £ 
generates density and entropy changes. However, this 
only accounts for five out of six possible phase-space di- 
mensions. Equation (|80f therefore implies a displace- 
ment direction that generates neither Sp, nor Ss, but still 
produces a momentum. The Lagrangian density only de- 
pends on two thermodynamic variables, and hence dis- 
placements associated with equation (|80|) represents a 
so-called "cyclic coordinate" familiar from Lagrangian 
particle mechanics. Each cyclic coordinate implies the 
conservation of an associated conjugate momenta. We 
show here that because of this invariance, the Lagrangian 
framework manifestly conserves energy and potential 
vorticity. 

First we define the following two quantities 



p ou 



fi = V xU 



(81) 



For simple Newtonian dynamics, U = u and O repre- 
sents the traditional vorticity; in a rotating frame 1~2 also 
includes contributions from the background "planetary" 
vorticity. Then defining the following quantity, 



A=— + ft xu + V(U 

ot 



u 



(82) 



variation of the action over tp from equation (|80[) implies 



Vs • (V X A) = 0. 

Defining the generalized potential vorticity 



GPV = 



P 



(83) 



(84) 



and then combining equations ([5]) & ([5]) and some vec- 
tor calculus identities gives conservation of GPV along 
trajectories of the flow 



d „ 



Vs „ fldC 

— -Vx ~5r 

p \p Ou 



= 0. 



(85) 



Equation (|85[) (for any C) represents a significant addi- 
tional advantage of the Lagrangian approach, especially 
when we begin to apply approximations to the fully com- 
pressible dynamics. 

In the case where the flow contains other frozen-in dy- 
namical quantities — such as magnetism in Appendix [Dl 
or chemical concentration — then variation in the direc- 
tion of equation (|50|) produces explicit changes in the 
Lagrangian density and equation (|85[) no longer holds. 

3.4. Lagrangian Derivation of the Euler Equations 

With equation ([73|) , we derive here the FC momentum 
equation with the Lagrangian density, 



C 



u 



— - e(p,s) - <j> 



(86) 



where equation ([7]) gives e = e(p, s). Computing specifi- 
cally each relevant term 



dC 



pu, 



du 

dC \u\ 2 p 
~d~ P = ~~ e ~p~ ~ 
dC 

it=p t - 

as 

Substituting into equation (|73l) . 

d{ P u) vm 2 

V • (puu) + p 



(87) 



(89) 



dt 

V 2 p 



pTVs = 0. (90) 



The last two terms in equation (|90[) simplify because of 
the Second Law of Thermodynamics (equation (J7JI ) 

TVs = Ve- -^Vp, 

p 2 

and the FC momentum equation (1221) follows in compact 
form 

Du „ „ 

p— + Vp + P V^ =0. 

For the FC equations, the energy adopts the expected 
form 

Furthermore, the Legendre transform of C with respect 
to p reduces to the thermodynamic pressure 



V = C-p—= P 
op 



(91) 



At this point, using Hamilton's principle to derive the 
momentum equations merely reformulates a well-known 
equation set. The utility of equation ([73"]) lies in that we 
need not use the specific Lagrangian density C in equa- 
tion (|86| . When considering low Mach number flows we 
also possess the tools to produce constrained dynamics 
that manifestly filter acoustic modes. Moreover, any con- 
strained dynamics derived from equation (|73[) will auto- 
matically achieve proper energy balance; this also leads 
to correct energy budgeting even if the constraint de- 
pends explicitly on time. 

4. CONSTRAINED MODELS 

We now use the Euler-Lagrange equation (|75|) to derive 
a momentum equation for subsonic constrained flows. To 
obtain constrained dynamics, we first recall the hydro- 
static reference state 



Vp + p G Vq 



0. 



As a fluid parcel moves slowly through a background 
stratification its internal pressure must equilibrate 
rapidly to approximately that of its surroundings. If the 
parcel does not remain almost in pressure balance then it 
could isotropically compress or expand without moving 
vertically. Rapid compressive collapse and/or expansion 
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generate the fast oscillations characterizing acoustic mo- 
tions. We wish to remove these modes from the dynam- 
ics. In this section we do so, and we will see that a succes- 
sion of more restrictive assumptions leads to the gener- 
alized pseudo-incompressible, anclastic, and Boussinesq 
equations in turn. 

4.1. Lagrangian Derivation of the Generalized Pseudo- 
Incompressible Equations 

The first and least restrictive constraint is that pres- 
sure fluctuations equilibrate rapidly and are small. We 
therefore impose the following (velocity independent, or 
holonomic) constraint on the dynamics 



C =p(p,s) -p (x) = 0. 



(92) 



Enforcing pressure equilibrium implicitly requires the 
fluid to move in some reduced subset of the phase space. 

The example of a roller coaster on a track provides 
an appropriate analogy to the fluid problem. The track 
(background pressure stratification) is strictly speaking a 
dynamical system in and of itself, and in particular sup- 
ports elastic flexural waves (acoustic waves). The roller 
coaster (subsonic swirling motion) couples to the very 
stiff (high sound speed) steel rails, which flex and push 
on the train as it moves along (fast acoustic restoration). 
But simply requiring the constraint (pressure equilib- 
rium) says that the steel in the track responds so rapidly, 
that effectively the train must slide along parallel to the 
rails. In such a system, the train's momentum is not con- 
served. The track must impart momentum to the train, 
otherwise the roller coaster would merely plunge straight 
downward (radially divergent pressure-driven collapse). 
If the initial energy of the roller car is sufficiently large 
compared to the potential energy in the rails (high Mach 
number situation), then the constraint no longer suffices, 
and the dynamics of the car and track progress in a 
fashion where it becomes difficult to distinguish between 
the separate dynamics of the two components. Ignoring 
here the possibility of a catastrophe (high Mach number 
flows), in a frictionless system energy remains conserved 
and we may derive the equations of motion via the ac- 
tion principle paired with a Lagrange multiplier and con- 
straint. 

Here we derive equations of motion for low-Mach num- 
ber subsonic flows using constrained Lagrangian analysis. 
For a system under constraint, the Lagrangian density 
becomes 

£->£- XC(p,s), (93) 

where equation (|92[) gives the constraint C, and A repre- 
sents the Lagrange multiplier, an additional independent 
variable of the action. Thus the full Lagrangian density 
becomes 



C 



P 



— e(p,s) -cf>) - XC(p,s). 



(94) 



The appropriate new term in the Euler-Lagrange equa- 
tion (|73|) assumes that stationarity of the action with 
respect to A enforces the constraint, 

dC 

— = -C(p,s) = -p + po = 0. (95) 

Stationarity of the action with respect to £ yields the 
equations of motion as in equation (|73|) . 



The strictly additive constraint in equation (|94l) im- 
plies the portion of the dynamical equations deriving 
from the original (unconstrained) Lagrangian density 
remain unaltered. Therefore, the constrained Euler- 
Lagrange equations assume the general form 



Du „ 

dC 



pV<f> = 



-Vlr.Ari + AVp, 



where we use 



dC 
dp 



dC 

ds 



(96) 



Defining the rescaled Lagrange multiplier 



P = T 1P X, (97) 
the constrained momentum equation now becomes 



Du _ 

Pj^ + Vp + pV4> = -Vp 



P 
Tip 



Vp. 



(98) 



Comparing equation with equation (|2"2")> . we see that 
equation (|98[) now contains an additional source of mo- 
mentum that prevents the dynamics from leaving the 
constraint manifold. 

The constraint equation (|9"2")l allows us to replace freely 
p — > pq, leading to the following constrained low-Mach 
number equations for momentum, density, and energy: 



Du 

P Dt' 
Dp = 
Dt 

V • u 



(p - Po) V<^> = 
u • Vpg 

rx(po,p)po' 

u • Vlnpo 

Ti(po,p) 



-Vp- 



p Vlnpo 

riOo,p)' 



(99) 
(100) 
(101) 



We dub this equations (|M]) - (|101[) the Generalized 
Pseudo-Incompressible (GPI) equations. The scaled La- 
grange multiplier p enforces equation (|101[) . In equa- 
tion (|9"9"j) we use equations (fT5"j) & (|9"2"|) to replace Vp — > 
Vpo = — PoV0, but the generality of equation (|98|) docs 
not require hydrostatic balance. In the limit of an ideal 
gas equation of state with constant 7, equations (1991) 
(fTUTI) reduce to the PI equations (Jp}-([1)|). Thus, the 
GPI equations extend the PI equations to more general 
equations of state. 

Most previous derivations of anclastic or pseudo- 
incompressible models follow from an asymptotic expan- 
sion of the full equations under the assumption of par- 
ticular scalings of the leading-order terms. The asymp- 
totic proce dure applies to many sy st ems under strong 
constraint (|Julien k Knobloclj[2007n . ISpieeel fc VeToliisI 
(|196C0 pioneered this approach with the Boussines q ap- 
proximation for a compressible fluid, iGoughl (|1969[ ) does 
this for t he an el astic models di scuss ed in P art I, while 
iDurranl (|1989( ). lAchatz et ail (|2010j ). and Klein et al.l 
(|2010f ) do this for the pseudo-incompressible model. In 
all of these derivations, the authors implicitly assume 
\p — pa\ <C Po, which may or may not imply a similar 
relation for density. 
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Our approach makes the same assumptions as previous 
authors. Our method differs in that we enforce the pres- 
sure constraint identically; i.e., non-perturbativcly. The 
main disadvantage of our approach is that our constraint 
tells us nothing quantitative about when our assumptions 
breakdown. We must apply physical reasoning (as in the 
roller coster analogy) and solution monitoring to deter- 
mine when acoustic and swirling motions loose their dis- 
tinct identity and thus when our derived equations lose 
their validity. The main advantage of our derivation is 
that it provides a general, systematic, framework for ex- 
tending the PI equations to include additional physics. 
Hamilton's principle gives a systematic method for gen- 
eralizing the PI equations to an arbitrary equation of 
state, while correcting the energy violations that emerge 
in the LM equations. 

We can apply the same Lagrangian framework to other 
important problems in astrophysical and geophysical 
fluid mechanics. For example, we derive an analogous 
model in Appendix [Dl that incorporates magnetism cor- 
rectly in a general framework. We note that any short- 
comings of the GPI equations derived here also apply to 
the previous sets of sound-proof equations discussed in 
2] In particular, the GPI equations reduce to the PI 
equations for an ideal gas and to the anelastic equations 
for an adiabatic background stratification. Furthermore, 
even if a filtered model gives lacklustre quantitative re- 
sults, the equations may still give excellent qualitative 
insight by eliminating complicating details. Finally, well- 
motivated subsonic models may guide the construction of 
numerical methods that take fast dynamics into account 
in an optimal fashion. 

The energy density in the GPI equations retains the 
same form as for the the FC equations, i.e., 



(102) 



where implicitly p = pq. 

However, we score more physical interpretation when 
computing the generalized pressure, 



C-p[ — e-4>- 



P _\(^P 
P 



dp 



(103) 



Simplifying equation (|103[) and using equation (|97[) we 
find, 

V=p + p = po+p (104) 

Comparing to equation equation (|104[) suggests 

that p indeed carries the closest possible correspondence 
to the pressure fluctuation p <-> p' = p — pg in the fully 
compressible dynamics. We show in the next section that 
the correspondence between the Lagrange multiplier and 
pressure perturbations arises naturally from considering 
enthalpy in the Lagrangian framework. 

For the GPI equations, the potential vorticity be- 
comes, 

,= f&-^.(V Xtt ), (105) 
where c p , T, olt, and T\ all implicitly assume p = 



Pq. Given equation (|105j) . equations (|M)) - (|101|) imply 
dtq + u • Vg = directly from the underlying Lagrangian 
structure, i.e., equation (|85l) . 

4.2. Enthalpy 

In §3.41 we express the Lagrangian density in terms 
of the fluid's internal energy, and hence its density and 
entropy. Here, w e provide an alte r nativ e derivation in 
terms of enthalpy. iSalmon fc Smith! (|1994|) used the same 
transformation in the Hamiltonian derivation a nonhy- 
drostatic pressure-coordinate model. This requires us to 
introduce the pressure as a canonical variable, in addition 
to and requires stationarity of the action with respect 
to this additional degree of freedom. This approach con- 
tains a number of advantages in helping to clarify the 
meaning of the Lagrange multiplier, and also points the 
way to more sophisticated approximation schemes in the 
future. 

Therefore, we make the following replacement 



p 

e(p,s) -> h(p,s) 

P 



(106) 



in equation (|86[) to produce 



£(u,p,s,p)=p— ph(p, s) - p(f> + p. (107) 

Stationarity of the action with respect to p exactly en- 
forces the equation of state given in equation (|10p , 



dC 
dp 



= -P 



dh(p, s] 
dp 



1 = 0. 



(108) 



In this sense, the pressure acts as a nonlinear Lagrange 
multiplier; similar to a generalization of XC(p, s,x) in 

We make the pseudo-incompressible assumption 



p^po(x)+pi, with |pi| < |p |. 



(109) 



Expanding the Lagrangian density in equation (|107l) to 
linear order in p\ gives 



C 



-p(h(p ,s) + (j))-pi 



P 



p(po,s) 



1 



-Po- (HO) 



Now pi acts as the Lagrange multiplier enforcing the 
constraint 



1 



p(Po,s) 



dh 
dp 



P=P0 



1 

p 



(111) 



Equation (jllll) requires some interpretation. The den- 
sity on the right-hand side of equation (jllip (represent- 
ing mass per unit volume) remains froze n into the flow , 
and evolves according to equation ([5]). IDurra n (1989) 
used the notation p* , to distinguish this from the density 
in the equation of state. The density on the left-hand side 
of equation (|111[) represents the equation of state density 
derived from enthalpy. Equation (|lllj) implies that the 
dynamical mass density must instantaneously equilibrate 
to the potential density of a parcel adiabatically trans- 
ported to a reference pressure level, po- This is math- 
ematically equivalent to the constraint p(p,s) = po(x) 
imposed in £14.11 

Furthermore, substituting equation f| 1 10[) into the gen- 
eral Euler-Lagrangc equation ((73)) gives exactly the same 
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result as equation (|98[) only with p\ replacing p. This 
underlies the correspondence between the Lagrange mul- 
tiplier and the small dynamic pressure fluctuations. Un- 
derstanding this correspondence becomes necessary when 
one wants to import data from a sound-filtered numeri- 
cal model into a fully compressible simulation. The fully 
compressible simulation will generate significant acoustic 
transients if the initial conditions start out of pressure 
balance. Using equation (|98[) with the extra momentum 
source relative to the LM equation (|47| , allows the pres- 
sure field to donate just enough momentum to the flow 
to keep it from oscillating rapidly if it were allowed. 

In principle, one could choose to carry the expansion of 
enthalpy around the background pressure to higher-order 
in p\. This would reintroduce linear acoustic dynamics 
back into the equations of motion. This gives an interme- 
diate approach between the FC equations and the GPI 
equations, and in principle, would allow a scale- by- scale 
filtering of acoustic modes according to some criterion 
such as their timescale. We leave these types of consid- 
erations for future work. 

4.3. Lagrangian Derivation of the Anelastic Equations 

We next tie our current results back to Part I, showing 
how the LBR anelastic equations naturally follow from 
an additional linearization around the background en- 
tropy profile. 

We point out that the constraint in equation (|110[) al- 
most constrains the density to give the background value, 
except for the entropy dependence. Therefore, if we make 
the replacement in the constraint term, 



Pi 



1 



Pi 



p(Po,s ) 



1 



(112) 



then stationarity with respect to p\ would yield p = po 
and V • {pou) = 0. Making the replacement in equa- 
tion (|112|) produces a version of the LBR anelastic equa- 
tions with a nonlinear buoyancy term deriving from en- 
thalpy term in equation (|110j) . However, we cannot jus- 
tify equation (|112p without admitting that 

p ps po s ps s (113) 

in some fashion. In particular, we must guarantee that 
equation (|113[) is consist with the entropy equation. 

In Lagrangian coordinates, the dynamics conserves the 
specific entropy of a fluid parcel along a given trajectory. 
That is, 

s(t,x(t,a)) = s (a), (114) 

where a = x(t = 0, a) represents the initial position of a 
fluid parcel. Therefore, defining an entropy perturbation 

s'(t, x(t, a)) = s (a) - s (x(t, a)), (115) 

the mean- value theorem implies that 

\s'\ < L 2 max|Vs (x)|, (116) 

X 

if we assume the background profile only depends on the 
height within of the layer. Here L z represents the total 
layer thickness. 

Therefore, if we assume a small background entropy 
gradient, then we may further linearize the Lagrangian 



density in equation (|110[) and obtain 

12 



£ ~ Po—^ PoToS-pi 



- 1 



Po 



Co 



(117) 



where £ = — Po{ e o — T s + 4>) gives a dynamically ir- 
relevant offset. Equation (|117p implies that we neglect 
energy sources from nonlinear terms in the thermody- 
namic variables. Extremizing equation (|117[) produces 
the compact form of the LBR formulation of the anelas- 
tic equations, 



Du 
~Dt 

Ds' 
Dt 

V • (p u 



where often 



Vzu = -s'VTq 
u • Vs = 

= o, 
Pi 

Po 



(118) 

(119) 
(120) 

(121) 



denotes a kinematic pressure variable, and p\ = P\jc\ 
gives the physical density perturbations associated with 
the evolution in equations (| 1 18|) (|12 1[) . This derivation 
remains valid only if s' remains small. By equation (|11GI) . 
this can occur in general only if Vso is small compared 
to ij 1 . This holds clearly in a nearly adiabatic atmo- 
sphere, where Vso is tiny. Alternatively, the assumption 
of small s' may hold (weakly) for motions with small 
vertical displacements compared to L z . 
For an ideal gas, 



VT = -^=* 



(122) 



gives the temperature gradient in the more common form 
for an adiabatic background atmosphere. This then puts 
the momentum equation in a more familiar form 



Du 
Hi 



Vzu = - 



s'9 



For a radiation-pressurc-dominatcd gas, c p 
ever 



VT = 



so 



(123) 



oo, how- 



(124) 



for the nearly constant background entropy so- As for 
the LBR equations of Part I, the LBR Equations (|118[) - 
(II 201) also conserve the energy 



E = po 



2 



s'Tr 



(125) 



and an associated free energy. Whilst the buoyancy 
term in equation (|118p may appear surprising on first 
glance, temperature represents the only physically mean- 
ingful quantity for nearly adiabatic stratifications that 
produces an energy when multiplied with entropy fluctu- 
ations. 

When compared to Part I, equation (|118[) implies a se- 
vere contradiction when attempting to model an isother- 
mal atmosphere with a "nearly adiabatic" background. 
This abuse of approximation becomes clear when one no- 
tices the complete vanishing of all buoyancy forcing in 
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equation Q118P with constant Tq. This implies that one 
may derive the isothermal LBR equations from a La- 
grangian that replaces To —4>/c p in equation (|117p . 
but this replacement would not follow from any sys- 
tematic approximation of the fully compressible dynam- 
ics, and would only apply for an ideal gas. Neverthe- 
less, e quat io ns (|M|) - ([TUl]) reduce to the anelastic equa- 
tions (|1 18|) (|120[) for a nearly adiabatic background. 

4.4. Boussinesq Equations 

At the outset, we only require a Lagrangian density 
pertaining to a fluid media, which depends on p and s, 
but not directly on the Lagrangian displacement as in 
the case of a solid. Other than the fluid requirement, we 
say nothing more specific about the equation of state — 
in particular we did not distinguish between liquids or 
gasses. The difference between liquids and gasses mainly 
derives from the strong independence of density as a 
function of pressure at fixed entropy, or T\ 3> 1. 

For example, in case of water under laboratory condi- 
tions, p ~ 10 5 N/m 2 , p w 10 3 kg/m 3 , c « 1500m/s, 
and therefore Ti 2.25 x 10 4 . H p = c%jg 225 km, 
However, over the temperature range of liquid water, the 
density can very by up to 10%, with a density m aximum 
at T 4°C allowing for penetrative convection (j Veronisl 
|1963[) . Therefore, nontrivial density backgrounds may 
exist in modelling liquids, but such systems may still fil- 
ter acoustic modes and may not undergo any significant 
volume changes. 

We therefore note that the GPI equations directly im- 
ply generalized Boussinesq equations in the limit 

Ti -> oo. (126) 



Dropping all terms involving Ti in equations ([99 |) ~ (|101[) 

gives 



Du 



r Dt 
Dp 

Dt 
V • u = 0, 



(p - Po) Vij 



-Vp, 



0. 



(127) 

(128) 
(129) 



Unlike the standard Boussinesq model, po may depend 
nontrivially on z, and large density fluctuation may oc- 
cur, especially in the presence of combined compositional 
and thermal effects. 

5. LINEAR WAVES IN ISOTHERMAL 
ATMOSPHERES 

As in Part I, we illustrate the properties of various sys- 
tems of equations with the dynamics of gravity waves in 
bounded plane-parallel atmospheres assuming an isothcr- 
mally stratified ideal gas. These simple atmospheres 
possess constant sound speed, buoyancy frequency, and 
density scale height. This makes computing eigenfre- 
quencies and eigenmodes for linear gravity and acoustic 
waves (when present) analytically tractable. This helps 
elucidate the differences between the various soundproof 
equations. These analytic results set the stage for the nu- 
merical experiments of |JSJ We use an ideal gas equation 
of state 



P= (cp - c v )pT = (7 - l)pe 



(130) 



with specific internal energy e, and 7 = c p /c v = 5/3. 

As in Part I, we define the velocity in terms of the 
vector displacement, 

u = -—, 

dt 7 



(131) 



which allows the simple integration of the linear density 
and pressure equations. 
We assume wavelike perturbations 



£ oc S(z) exp (iu)t — imx), 



(132) 



where x represents the horizontal coordinate, and m its 
associated wavenumber, and 'B(z) gives the vertical struc- 
ture depending on the specific model. 

In a hydrostatically balanced (equation (fT5"jl ) isother- 
mal atmosphere 



Vlnp = Vln/3 = - — 
tl 



19 z 



(133) 



where H is the (constant) pressure or density scale 
height, and 

JPo 
Po 



(134) 



is the (constant) sound speed. The Brunt- Vaisala fre- 
quency N is 



N 2 



g ■ Vfo = (7 - 1) .9 

H 



(135) 



c P 7 
where Vso is the background entropy gradient. 

5.1. Fully compressible waves 

The well-known solution for the full y compressible 
equations exists in several textbooks (e.g., Lighthill 1978) 
and in Part I. In a bounded or unbounded atmosphere, 
vertical eigenfunctions of the take the following form 



cxp Oh) 



sin(fcz + 8), 



(136) 



where k and <5 give an arbitrary vertical wavenumber and 
phase. For finite domains with impenetrable boundaries, 
only quantized modes 



7rn 
k = —, n 

Liz 



1,2, 



(5 = 0, 



(137) 



ensure that vertical motions cease (£ z = 0) at z = 0,L Z . 
Equation (|136|) implies that the frequencies of gravity 
and acoustic modes remain purely real, i.e., 



1 



AH 2 



(138) 



The quadratic nature of equation (|138[) in u> 2 provides 
the two distinct acoustic and gravity wave branches. We 
write the full solution to both branches of equation t|138[) 
in the form, 



1±4 1- 



4w 2 



where 



m 



AH 2 



(139) 



(140) 
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and 



k 2 + m 2 + jw- 



-JVrf 



(141) 



Equations (|140j) & (| 141[) represent simple approxima- 
tions to the behavior of equations (|138[) & (|139|) in the 
high- and low-frequency limits respectively. 

Importantly, the linear waves conserve the leading- 
order quadratic energy in equations (|30l) - (|32[) . 

5.2. Pseudo-incompressible gravity waves 

We begin our linear analysis of the sound-filtered sys- 
tems with the energy-conserving PI equations. As we 
show in the GPI equations (j9"9"j) - (fTU01) reduce to 

the PI equations (j3"5)) - (|4T))) for Ti = 7. Linearizing the 
pressure, momentum and density equations around the 
background state gives 



p uj 2 $- 



-PoV 



pig, 



P\/po = -£ ■ Vlnpo - V • £ 
V •£ = -£• Vln/3 , 



(142) 

(143) 
(144) 



where /3q = p ^/ 1 ( equation (|20|) ). We combine the den- 
sity equation (|143[) and pressure equation (|144|) to give 

Po 



Po 



(145) 



Thus the buoyancy term in the pseudo-incompressible 
equations remains identical to the buoyancy term in the 
anclastic equations considered in Part I. 

The divergence of the horizontal momentum equation, 
and the pressure constraint determines the pressure fluc- 
tuations in terms of the vertical displacement, 



VV = p uj' z Vj_ • £j 



u 2 Po d(Mz) 
' Po 



dz 



(146) 



Combining the vertical momentum equation f| 142[) and 
the buoyancy equation (|145[) for linearized waves, we ob- 
tain 



2 V 2, Po d 

Po oz 



Po 



(147) 



where we use equation (|135|) . Collapsing the entire sys- 
tem into a single second-order equation for vertical dis- 
placement gives 



where, 



2 /?Q_9 

Po 



« = N 2 Vl£ z , 
Po d(Mz) 



PI dz 



(148) 



(149) 



is a negative-definite self-adjoint operator with respect 
to the integration weight function po(z). Sclf-adjointncss 
follows from the identity 



Po V £>pi£ dz = - 



Po d(P oV ) d(M) 



(150) 



assuming r] 
domain. 



Pq dz dz 
on the endpoints of the integration 



For the isothermal background in H5.ll the same ver- 
tical eigenfunction profile in equation (|136[) provides the 
real- valued dispersion relation 



2 



2 M ,2 , (7-2)' 



4 7 2#2 



m 2 NL 



(151) 



which closely resembles equation (|141j) . but with a differ- 
ent large-scale cutoff. We note that the correspondence 
between PI and FC eigenfunctions breaks down in spheri- 
cal geometries. In spherical geometries, the radial struc- 
ture of the FC eigenfunctions depend explicitly on the 
spherical harmonic degree, I, but not in the PI approx- 
imation. However, many of the fundamental qualitative 
properties developed here remain even in those systems. 

In a monatomic ideal gas with 7 = 5/3, equation (|151|) 
implies a large-scale cutoff frequency 1/(107J); the cutoff 
vanishes when 7 = 2, and approaches 1 / (2H) for larger 7. 
The vanishing of the large-scale cutoff for 7 = 2 actually 
corresponds to similar behavior in the fully compressible 
dispersion relationship for any m as k — > 0. More gener- 
ally, for any background, equation (|15ip gives an optimal 
approximation to the FC frequencies in the sense that 



wpi| = O (m 



(152) 



More gencrically, the error in the difference between the 
approximate and fully compressible frequencies scales as 
~ m . 

Figure [1] shows the dispersion relationship for the PI 
equation (|15ip . in the limit kH — > 0, and compares 
this to the equivalent modes in FC, LM, and all other 
anelastic models considered in Part I. Figure Q] clearly 
shows the behavior in equation (|152[) relative to other 
models for large mH. At low horizontal wave number 
m, the reduced large-sale cutoff in the PI dispersion re- 
lationship (equation (|151jl ) leads to higher frequencies 
than those obtained from the FC equation (|139|) . For 
particular comparison, we also plot the dispersion rela- 
tionship of the energy-conserving LBR anelastic equa- 
tions from Part I. Generally, frequencies obtained from 
the LBR equations match the Euler frequencies well at 
low frequencies but approach uj ~ Nq slowly for large 
m. The differences in frequencies obtained in all energy- 
conserving equation sets become smaller as kH increases. 
We also show all of the equation sets considered in this 
paper and in Part I in Figure [T] for the same isothermal 
atmosphere. The PI, FC and LBR equations conserve 
energy, while the LM and ANS equations do not. In 
general, the RG equations do not conserve energy, ex- 
cept in the special case of an isothermal atmosphere (see 
Part I). At large mH, the ANS, LBR, LM, PI and FC 
equations converge to the Brunt- Vaisala frequency Nq, 
while the RG equations are too large by a factor of -Jj. 
Table [1] gives the dispersion relationships for all equation 
sets in both papers. 

Figure [2h shows the first five vertical modes fci-fes of 
the dispersion relationship for the PI, LM and FC equa- 
tions, where equation (|137j) gives k with L z = bH. Simi- 
lar to Figure [H the frequencies of gravity waves in the PI 
equations remain somewhat larger than the FC equations 
at low horizontal wave number, with the discrepancy 
most notable for long- vertical- wavelength modes (e.g., 
fci ) . Figure [2j> shows the vertical eigenfunction struc- 
ture for the ki mode. This mode, and the stratification 
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Figure 1. Dispersion relationships for waves with kH — > for all 
equation sets considered here and in Part I (and see Tabled). Here 
we also show the sound wave branch (black, labelled "acoustic") of 
the exact solution to the full Eulcr equations. 



of this atmosphere, correspond to the numerical results 
we present in ^Bl 

As in Part I, we normalize these eigenfunctions with 
an amplitude A, such that 



A 



r 5H 

Jo e ~ 



Z / H dz 



2 1 



-5e ' 



(153) 



In the PI and FC equations, e = 1, and Api w 1.59, for 
7 = 5/3 whereas in the LM equations e = (7 + l)/7, and 
A hM « 2.00. 

5.3. Low Mach number gravity waves 

Finding linear cigenfrequencics in the LM equa- 
tions (j47 |) -(|49 | amounts to the same procedure as in 
Now however the linearized momentum equation is 



- p a u) 2 £, = -Vp' + pig. 



(154) 



The fundamental differences in the pressure term 
between equations (|142[) & (|154[) lead to the non- 
conservation of energy in the LM equations and instead 
to the conservation of pseudo-energy for linear waves. 

Using the linearized equations (|143p & (|144[) for den- 
sity and pressure, and equation (|1 54[) for momentum, we 
collapse the LM system into a form equivalent to equa- 
tion ([148]) . 



LO 2 (V? 



2>lm) 



where, 



V 



LMtz " p dz 



Po dz 



(155) 



(156) 



is a negative-definite self-adjoint operator with respect 
to the modified integration weight function, 



Po(z) = po(z)fio(z). 



(157) 



That is, assuming r\ = £ = on the endpoints of the 
integration domain, 



PoVVtutdz 



p d(p ri) d(M) 




Figure 2. Properties of waves in bounded atmospheres, (a) 
Dispersion relationships for PI equations (blue, solid), LM equa- 
tions (purple, dashed) and full Euler FC equations (black, solid) 
for waves with vertical wavenumbers k±-k^, as labelled, plotted 
against scaled horizontal wavenumber mH. The numerical simu- 
lations of [[6] use &2 and mH fa 1.26. (b) Vertical eigenfunc tions 
for the PI and LM equations for the &2 mode with equation 1 11531 
giving the normalization. Eigenfunctions of the FC equations are 
identical to the PI equations and are not shown. Though the 
frequencies of the LM equations arc reasonably similar to those 
obtained in the energy conserving PI and FC equations, the LM 
eigenfunctions differ significantly. 

For an adiabatic background atmosphere, po oc p^ and 
the LM equations do not reduce to an anelastic system. 
This implies that the LM equations would not give cor- 
rect growth rates and vertical structure for convection in 
an unstably stratified system. 

For an ideal gas, equations (|155[) - (|157[) require the ver- 
tical eigenfunction 

E(z)=cxp(^^-^jsm(kz + 5), (159) 

to guarantee real- valued wave frequencies, 

(7 -I) 2 



2 

W LM 



mf + k + 



4 7 2#2 



m 2 N^. (160) 



As with the ANS equations in Part I, a serious problem 
lurks with this choice, since the kinetic energy density 
scales as 



Pou oc exp 



1 h)' 



(161) 



PI 



dz dz 



dz. (158) Therefore, the kinetic energy of the waves grows ex- 
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Table 1 

dispersion relationship comparison 



System 


(u/rnNo)- 2 = 


eq 


FC 


see text 


fL39l) 


PI 




JTsTt 


LM 




<fT60l 


ANS 




(PI: 56) 


LBR 




(PI: 62) 


RG 




(PI: 67) 


Note. — We quote the ANS and LBR 
dispersion relationships from Part I (PI). 



poncntially with height. From the FC equations, we 
know that this kinetic energy density should remain con- 
stant with height in an isothermal atmosphere. This 
discrepancy arises from the fact that the LM equations 
fail to properly conserve energy, but rather conserve the 
pseudo-energy given in equation (l63l) . 

Figure [2]a shows the dispersion relationship with real 
uj given in equation (|160[) for an isothermal atmosphere 
with L z — 5H. Like the PI equations, the LM equa- 
tions give slightly higher frequencies than the fully com- 
pressible equations. A much larger problem however 
rests with the eigenfunctions of the waves. Figure [2^ 
shows the eigenfunctions in an isothermal atmosphere 
with L z = 5H. Equation (|153[) provides the normaliza- 
tion, and gives the response of waves in the LM equa- 
tions generated from the same initial density perturba- 
tion used for the PI equations. Clearly, the eigenfunc- 
tions of waves in the LM equations acquire an excess 
amplitude compared to the correct eigenfunctions given 
in equation (|136j) . 

6. SIMULATIONS OF GRAVITY WAVES IN A 
BOUNDED ISOTHERMAL ATMOSPHERE 

In this section we demonstrate through direct numeri- 
cal simulation the properties of the PI and LM equations 
described in gSHEll &H We use the MAESTRO code 
to simulate the propagation of internal gravity waves in 
a bounded isothermal atmosphere, implementing both 
the PI and LM equations (see Appendix [E] for more de- 
tails). Although MAESTRO supports adaptive mesh 
refinement (AMR) and an evolvin g base state, we use 
neither of these ca pabili ties here ( Almgren et al.1 [20081 : 
iNonaka et aDIMol l20lll ). In our implementation of the 
PI equations , we take an ideal gas equation of state, 
assume that Ti(po,p) = 7 = 5/3, and solve the nonlin- 
ear PI equations ((35])-(|in|) and LM equations (|4T )) - (|4!?)) . 
In both cases we do not include any explicit diffusivi- 
ties, and rely instead on numerical diffusivitics within 
the MAESTRO code itself. 

We use an isothermal, hydrostatically balanced back- 
ground state with constant scale height, H , and gravita- 
tional acceleration, g, such that 

Po(z) = Prei exp(-z/H), p (z) = gHpo(z). (162) 

Although the PI and LM equations allow arbitrarily large 



density perturbations, both assume small pressure per- 
turbations, and thus Po(z) cx exp(—z/^fH) remains fixed 
throughout our simulations. In this paper we focus on 
a set of 2D simulations in a Cartesian domain with size 
L x x L z , where L x = L z = 5H. 

For any given background variable, say p, defining the 
number of scale heights of this variable proves helpful in 
discussing the degree of stratification in our simulations, 



In 



Po|z=o 

Po\z=L M 



(163) 



Therefore, n p = n p = 5 in this isothermal atmosphere. 
By comparison, for the radiative zone of the Sun, n p sa 
6.6. 

In the LM equations, /3q scales the pseudo-density ac- 
cording to equation (|157p . The number of pseudo-density 
scale heights rip is 



Up 



1 

n p + -rir 

7 



(164) 



where np and n p represent the number of /3 and pres- 
sure scale heights respectively. When n p = 5, np = 3, 
and rip = 8. Thus, the gravity waves in the LM simu- 
lations experience the effects of stratification much more 
strongly than the gravity waves in the PI simulations. 

We conduct the simulations in a Cartesian box with 
resolution (N X ,N Z ) = (512, 512). The box has a periodic 
horizontal direction, with impenetrable vertical bound- 
ary conditions at the top and bottom. Since we do not 
include explicit viscosity or thermal diffusion, we do not 
require any further boundary condition at the lower or 
upper boundaries. We initialize the various simulations 
with one of three small density perturbations, 



/2ttx\ . (2nz\ 
dppi = A a sin I I sin I — I exp 

(2-kx\ . (2-kz 



Splm = A sin I — — I sm I — - I exp 



-2H1 ' 
(7 + 1)* 



2jH 



„ . [2ttx\ . (2-kz 
OPGW = A sm I I sm I — 



(165) 
(166) 
(167) 



where Aq = 10 _6 p ro f, and p rc { gives the density at 
z = from equation (|162[) . The perturbations Sppi and 
SpLM correspond to the n = 2 eigenfunctions derived in 
0] for the PI and LM equations respectively, with an 
m = 10 horizontal perturbation. Importantly, each PI 
(LM) eigenfunction projects onto many LM (PI) eigen- 
functions, so the PI (LM) eigenfunction excites many 
modes of the LM (PI) equations. The third perturbation, 
Spow, projects onto many modes from both models, and 
thus excites a broad band of modes in either the PI or LM 
equations. We normalize the background density to one 
at the bottom of the domain, i.e., p re f(0) = 1. Because 
the density perturbations produce small velocity ampli- 
tudes, the default velocity-based CFL time-step condi- 
tion insufficiently captures the Brunt- Vaisala timescale 
of the gravity waves, where tbv = 

Nq 1 . Wc thus re- 
strict our time steps such that At < 0.016A r " 1 , which 
accurately resolves the relevant waves. 

The analysis in fJS] shows that the vertical eigenfunc- 
tions differ for the two sets of equations. In Figure [3] 
we plot the rms value of the density perturbation at two 
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(a) (b) 
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Figure 3. The rms density perturbation at time t = (blue) and 
time t = 600Nq 1 (red). In (a, b) the simulations employ the PI 
equations, and in (e, d) the simulations employ the LM equations. 
Cases (a, c) each use the initial perturbation i5ppi, (b, d) use <5pLM- 
Cases (a) and (d), show no change in form with time, indicat- 
ing the perturbations coincide with eigenfunction of the respective 
equations. In (6) and (c) the initial perturbation differ from eigen- 
functions of the relevant equations and so instead excite a band of 
modes, which interfere at later times. 

times, t = and t = WONq 1 , in both (PI and LM) 
types of simulations using initial density perturbations 
based on each eigenfunction in equations (|165[) & (|166f) . 
We compute the rms perturbation by taking the square 
root of the horizontally averaged potential energy density 
(equation (|60])) times 2JVq /(poff 2 )- In each set of equa- 
tions, when we perturb the background with the cor- 
rect eigenfunction, the perturbation remains the same 
after an integer number of periods. We demonstrate 
this for the PI equations with initial perturbation dppi 
(Fig- Eta)), and the LM equations with initial perturba- 
tion 5/5lm (Fig. E[d))- Otherwise, the perturbation ex- 
cites a broad band of modes which interfere, and the rms 
amplitude changes in time (Fig. [3^0, c)). This confirms 
the predicted differences in the eigenfunctions of the PI 
and LM equations, as given in equations (|1 36[) & (|159[) , 

Another significant difference between the two sets 
of equations is that the PI equations conserve energy, 
whereas the LM equations do not. Although the LM 
equations violate energy conservation, they still contain 
a pseudo-energy quadratic invariant in the linear regime 
(equation (|63|) ) which differs from the true energy by a 
factor of /3n (see discussion in §2.4j) . 

We plot temporal traces of volume-averaged energy 
and pseudo-energy (equations (|58]l & (|63jl) for simula- 
tions implementing the PI equations and the LM equa- 
tions in Figure |H We integrate the equations until 
t = 1000-AT^ , and use <5pcw a s an initial perturbation, 
which excites a broad band of modes for both simula- 
tions. Figure HJa demonstrates that the PI equations 
conserve energy, whereas the LM equations clearly vio- 
late energy conservation (Figure Hb,e). Instead, the LM 
equations conserve the pseudo-energy (FigureUJi), which 
the PI equations do not conserve (Figure H^,/). Because 
both simulations use the the same initial perturbation, 
3pGW, the two simulations initially contain the same en- 
ergy, and the same pseudo-energy (Figure HJs,/). 

Figure His,/ show the long-time variation in energy and 
pseudo-energy. To quantify the deviations in these quan- 
tities, we define the energy and pseudo-energy variation 



as 

AE=(6E)/(E), (168) 
APE=(SPE)/{PE), (169) 

where (•) denotes a temporal mean and 6 denotes the 
standard deviation in time for the entire dataset from < 
t < 1000 Aq~ l . In the PI simulation depicted in Figure [4] 
AE w 0.035%, whereas APE 32%. Conversely, for the 
simulation implementing the LM equations, AE f=a 36% 
and APE 0.035%. 

Slowly accumulating errors within the explicit time- 
stepping method of MAESTRO leads to the slight devi- 
ations from perfect energy conservation in the PI simu- 
lation and perfect pseudo-energy conservation in the LM 
simulation. The small growth in the "conserved" quan- 
tities becomes even smaller when we decrease the size of 
our maximum time step. 

We compute a variety of simulations using the LM 
equations with simulations domains of various sizes, and 
with different initial perturbations. As in Part I, we find 
that the degree of energy non-conservation, measured us- 
ing AE, scales about linearly with n p . That is, larger 
n p leads to larger the variations in E. In the simula- 
tions shown in Figure 21 the initial perturbation contains 
only one horizontal mode. Including multiple horizontal 
modes lowers AE, since modes with different frequencies 
combine incoherently This does not mitigate the lack 
of local energy balance. In all cases, the pseudo-energy- 
conserving LM waves tend to have larger amplitudes than 
the energy-conserving PI waves. 

7. SUMMARY AND CONCLUSIONS 

In the last few decades, the fluid dynamics and mag- 
netohydrodynamics of stellar and planetary interiors has 
become amenable to direct numerical simulation. Al- 
though the simulations do not capture all the realistic 
parameter ranges — they tend to overemphasize the role 
of diffusive processes, for example — many approximate 
treatments such as the mixing length theory of convec- 
tion, Eddington-Sweet circulation, turbulent viscosity, 
and mean-field dynamo theory can now use numerical 
solutions of the fluid equations to calibrate useful param- 
eterizations. These advances make it possible to probe 
many dynamical phenomena, including penetrative con- 
vection, chemical mixing, angular momentum transport, 
and stellar and planetary dynamos, with unprecedented 
realism. 

Many interesting problems involve the coupling be- 
tween radiative and convection zones. Penetrative con- 
vection gives one example; the degree of overshoot affects 
dynamo-related phenomena such as magnetic pumping, 
chemical evolution phenomena such as depletion of light 
elements, and dynamical phenomena such as the gen- 
eration of gravity waves. Angular momentum balance 
generally also tends to couple radiative and convection 
zones; for example, red giant cores should spin up as they 
contract while red giant envelopes should spin down, but 
can angular momentum flow between these two zones? 

Because acoustic timescales in stars are typically very 
short compared to other dynamical timescales, direct nu- 
merical simulation of stellar interior dynamics is only 
practical if acoustic waves can be filtered out. Similar 
issues arise in the study of planetary atmospheres and ac- 
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Figure 4. Temporal evolution of energies and pseudo-energies in simulations solving the PI and LM equations with <5pg w as an initial 
perturbation, (a) Energy and (ft) pseu do-energy in the simulation solving the PI equations. The total energy (E, equation (1 58 1) , in black) 
and pseudo-energy (PE, equation J63l) arc divided by two to high ligh t the fluctuations bet wee n kinetic energy K (equation II59I I. in red) 
and potential energy U (equation (160 6 . in blue) or PK (equation | |64| |) and PU (equation H65H ) respectively. The PI equations correctly 
conserve total energy, but not the total pseudo-energy, (c) Energy and (d) pseudo-energy in the simulation solving the LM equations. The 
LM equations do not conserve the total energy, and instead conserve the total pseudo-energy. Wc plot the total energy (e) divided by two 
and total pseudo-energy (/) divided by two for the LM simulation (solid line) and the PI simulation (dashed line), over a longer timescale. 
The vertical dotted lines show the time range (SbONg 1 , 950NQ 1 ) depicted in (a-d). Avera ging over IOOCWq , we have that AE 0.035% 
for the PI simulation and APE fs 32%. Over this same timescale, AE « 36% for the LM simulation, but APE pa 0.035%. 
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cretion disks. Over the years, many ingenious and beau- 
tiful treatments have been developed to work around this 
problem. Two such approximations, the anelastic and 
pseudo-incompressible models, are valid for adiabatically 
stratified systems, and thus should be quite accurate in 
stellar convection zones. However, in recent years, inter- 
est in modeling stellar dynamos and stellar internal rota- 
tion, circulation, and mixing have also demonstrated the 
need for models which are valid in stably stratified, radia- 
tive zones. Additionally, modeling double-diffusive buoy- 
ancy such as fingering convection and semi-convection 
requires capturing both significant turbulent and gravity 
wave dynamics. 

In Part I we studied the anelastic equations and showed 
that one particular version (the LBR equations) was par- 
ticularly versatile because it conserves energy. In this 
paper we extend our treatment to two consider sets of 
pseudo- incompressible models (PI and LM). 

Since there are many methods for achieving acoustic 
filtering, which one works best? It has long been known 
that the fully compressible Euler equations can be de- 
rived from Lagrangian analysis. In <2] wc showed that 
approximations such as pseudo-incompressibility can be 
implemented by a Lagrangian analysis with the approxi- 
mation scheme enforced through constraints. The result- 
ing equations conserve energy in their ideal form and arc 
thus preferable over nonconservative formulations. As an 
important illustration of the technique, we extended the 
pseudo-incompressible approximation to systems with a 
generalized equation of state. This extension will enable 
treatment of gas-radiation mixtures, partially degenerate 
systems, and other thermodynamically complex states 
frequently encountered in stellar astrophysics. We call 
this model the generalized pseudo-incompressible equa- 
tions (GPI equations; JOJ. More generally, we pro- 
vide a systematic framework for producing approximate 
equations that filter fast dynamics. In Appendix [D] wc 
apply this to the magnetohydrodynamic version of the 
GPI equations. In the future, we expect to apply the 
Lagrangian framework to many other similar problems 
in astrophysical fluid dynamics such as rapidly rotating, 
strongly sheared, and/or multiple-component systems. 

Simple wave propagation problems provide a good 
test of the properties of different equation sets. In §5\ 
we derive and illustrated the behavior of gravity waves 
in the fully compressible (FC), pseudo- incompressible 
(PI/GPI), and low Mach number (LM) equations. Wc 
demonstrate that in an isothermal atmosphere, the en- 
ergy flux in FC and PI/GPI gravity waves remains con- 
stant with height, but that the LM gravity wave flux di- 
verges with height. These properties follow from energy 
conservation and non-conservation, respectively. 

In SJS] we conduct direct numerical simulations of lin- 
ear gravity waves in isothermal atmospheres using two 
versions of the MAESTRO code; the standard version, 
which employs the LM equations, and a modified ver- 
sion which solves the PI/GPI equations. We showed in 
Figures 4 and 5 that the cigenfunctions, energy conserva- 
tion properties, and pseudo-energy conservation proper- 
ties are as expected. Thus, we show that implementation 
of the PI/GPI equations is achievable and we hope that 
it will prove fruitful. 

In this paper and in Part I we worked only with ideal 
equations which omit viscosity, diffusion, and any other 



non-adiabatic heat transport effects. While it might be 
argued that energy conservation is only important in 
ideal systems, we note that if an ideal system does not 
conserve energy, there is no guarantee that non-ideal ef- 
fects will decrease energy; they might cause spurious in- 
stabilities or other poorly controlled behavior. Thus, we 
recommend always working with equations which con- 
serve energy in the ideal limit. 
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APPENDIX 

A. RADIATION-GAS EQUATIONS OF STATE 

Astrophysical fluids can depart significantly from the 
ideal gas equation of state assumption. In this appendix 
we show how to implement a more complicated equa- 
tion of state into the pressure-constrained GPI equations. 
That is, we show how to compute Ti(po(x), p) in a non- 
trivial situation. 

As an example, we present a mixed radiation and gas 
model, which commonly applies in astrophysical systems 
ranging from the interiors of massive stars to accretion 
disks around black holes. In this case, 

Pg = —pT, p r = (7r - l)aT^ T , p = p g + p r (Al) 

7Tip 

specify the (ideal) gas, (blackbody) radiation, and com- 
bined total pressures respectively (/cb, mp, and a repre- 
sent the Boltzmann constant, mass per gas particle, a 
the radiation constant respectively). For the radiation 
adiabatic exponent 7 r = 1 + l/n where n is the dimen- 
sionality of the space (e.g., j r = 3/2 in 2D and j r =4/3 
in 3D). 
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Equations for the internal energy e and entropy s follow 
from equation ([?])■ 



pe 



Pr 



Vg 



7r ~ 1 Ig - 1 ' 



log- 



(A2) 
(A3) 



The temperature parameterizes the various thermody- 
namic variables, implicitly giving p = p{p,s). In this 
case the first adiabatic index mixes contributions from 
gas and radiation, i.e., 

Ti = 

7g(7r - l) 2 Pg + lr{lg ~ l)Pr ((-for ~ l)Pg + 7rPr)) 
(Pr +P S ) ((7r - 1) 2 P 3 + 7r(7s - l)Pr) 

32 - 24A - 3A 2 



24- 21A 



(A4) 



where X = p g /(p g + p r ), l g = 5/3, and 7 r = 4/3 for 
the last expression in equation (|A4|) . E quation (IA4I) 
for 7,. = 4/3 matches the result given in | Cox fc Giulil 
(|1969f ). which also presents a wide range of astrophys- 
ically relevant equations of state. When gas pressure 
dominates Ti w 7 9 = 5/3 (for a monatomic gas in 
three dimensions). When radiation pressure dominates 
Ti « 7 r = 1 + 1/n = 4/3 (in three dimensions). Im- 
portantly, when neither gas nor radiation pressure dom- 
inate, then Ti depends on both {T, p}, or implicitly on 
any other independent pairs of thermodynamic variables 
such as {p, p}, or {p, s}. 

When applying the pressure constraint p = po, we must 
invert the implicit relationship 



, . fcn aT 
p (x) = — P T+ — 
m D 3 



(A5) 



by solving equation (|A5I) for T = T(po(x), p). The com- 
plete formula requires the roots of a 4th-ordcr polyno- 
mial, but one may easily numerically solve an equation 
of the form y = x + x 4 for x in terms of y. Finding T de- 
termines p g = p g {p(){x), p) andp r = p r (po(x), p) individ- 
ually. Substituting into equation (|A4[) gives Ti(p (x), p). 

B. FREE ENERGY 

In this appendix we derive the desired entropy function 
to produce a leading-order quadratic energy invariant. 
To start, we define 



F - 



p\u\ 



p\u\ 



+ p(e(p,s) + 0-f(s)) 
+ A(p, s) -p . 



where 



A(p, s) = po + p (e(p, s) - f(s) + . 



(Bl) 



(B2) 



represents the available potential energy. For notational 
convince, we define the perturbation variables 



si 



s , p 1 = p-p . 



(B3) 



Expanding the available potential energy to second order 
in the perturbations, 



A(p,s) = A 

d 2 A 



PlSl 



dpds 



dA 
dp 
s\ d 2 A 



si- 



dA 



p\ 3 2 A 



2 dp 2 



o 



2 8s 2 



of, sf, sxpx, pis\, 



(B4) 



where |o denotes evaluating at p = po and s = so- Re- 
quiring 

^(po,so) = e o + ^ + ^-/( So ) = (B5) 
op Po 

dA 

ir (po,s )=po(T -f(s )) = 0, (B6) 
us 

ensures that the linear (non-sign-definite) terms in equa- 
tion (|B4[) vanish, and the quadratic terms produce the 
first n on-t rivia l co ntributions to the free energy. Equa- 
tions (|B5j) & (|B6[) imply the following functional rela- 
tionships 

f(a {x)) = ho(x)+<f>(x), f'(s (x))=T (x), (B7) 

for background enthalpy, ho, and temperature, To. Equa- 
tion (|B7j) only allows a solution if the two relations do 
not contradict each other. Checking this, 



d(/io + 0-/(so))=d/i o 



dpo 
Po 



d/(s ) 



= (T -f(s ))ds = 0, (B8) 
the enthalpy form of equation ([9]) imp lies compatibility 



between the two relatio ns in equ ation (|B7[) 

Assuming equations (|B5j) & (|B6|) . we compute the 
leading-order term in the available energy, 

A(p ,s )=0. (B9) 

To compute the quadratic terms in equation (|B4|) ■ we 
first define the total background pressure differential. 



dpo = -pod(j) 



Op dp 

— (p ,s )dp + -^-{po,s )ds , (BIO) 

op OS 



and differentiation with respect to the gravitational po- 
tential function, 



Together, 

dp 



ds 

Therefore. 
d 2 A 



(Po,sq) = - 



dSr 



PO 



dp 



dp 2 
d 2 A 



(po,s ) 



dpds 
d 2 A 



{Po, so) 



i 

PO 

I dp 
= —Tr{Po,s ) 
Po os 



ds 2 



(po,s ) = Po /"(s ) 



dT 



(Po,s q ) 



Using 



/'(so) = n 



-g^(P0,S0), 



(Bll) 

(B12) 

(BI3) 
(B14) 
(BI5) 

(B16) 
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and differentiating with respect to 



dT 



<1sq 1 dpo dp 



fisa) ._ M ^ = ^ feM . (B17 ) 

To linear order, density and entropy perturbations gen- 
erate pressure perturbations via 

dp 

Pi = Pi Co + si — (po,s ) + 0pl,sl,p 1 s 1 . (B18) 

Putting everything together, we obtain the available po- 
tential energy to quadratic order in terms of pressure and 
entropy perturbations, 



A(p,s) 
where 



h = si 



ds ' 



(B20) 



The background stratification deter mines the sign of the 
second coefficient in equation (|B19j) . and therefore deter- 
mines the static stability to linear perturbations. 

C. BUOYANCY FREQUENCY AND ENTROPY 

Simply from the properties of differentials, the fol- 
lowing two relations hold for any three thermodynamic 
quantities, 



db I A da 



db J A dc J A da 



Computing the entropy differential 



(C2) 



equation (|C1[) implies 

_Udp 

p\ds 



1 / dp\ p ( dp\ f dp\ dp dp 

P\ds) p \dp) s p p 
dp dp 



Fi(p,p)p p 



(C3) 



Carrying this further, and using commutativity of second 
partial derivatives (Maxwell relations) 

1 (dp\ 1 (dv\ 1 d 2 h 1 (dT 



p\dsj p v \ ds J v dsdp v \dp 

(t~)(w~) = ~~(t) °° 

v\ds J p \dp J T vcp \ dp ) T 
T d 2 g T ( dv \ _ Ta T 



vc p dTdp vc p \ dT J c p 
Therefore, 

Ta T — = dP - — 
c p Tx(p,p)p p 

1 (dp 



where 



P \dT 



(C4) 



(C5) 



(C6) 



defines the thermal expansion coefficient at constant 
pressure, and 



(C7) 



defines the specific heat capacity at constant pressure. 
The Brunt- Valsala frequency in terms of the entropy gra- 
dient and gravity follows from equation (|C5|) . 



N * = -Ta T i^±=g.(*£-p). (C8) 



FipJ 



For an ideal gas, ar = 1/T, Ti = 7. 



D. CONSTRAINED 
MAGNETOHYDRODYNAMICS 

Many astrophysical fluids arc magnetized plasmas, 
and can be modelled by coupling the Euler equations 
with electrodynamics to derive the magnctohydrody- 
namic (MHD) equations. Assuming the fluid is a per- 
fect conductor, the magnetic field remains frozen into the 
fluid, and evolves according to the induction equation 



dB 

^— = V x (u x B). 

dt x ' 



(Dl) 



The magnetic field influences the velocity evolution 
through the Lorentz force, /i^ 1 (V X B) X B, which ap- 
pears in the momentum equation. 

Linear perturbations around a background flow and 
magnetic evolve according to 



dSB 
~dT 



V X (dux B + ux SB). 



(D2) 



Parameterizing perturbations using a displacement vec- 
tor ^ defined in equation (|70[) . one can integrate equa- 
tion (lD2l) such that 



6B = V X (£ X B). 



(D3) 



The special case wher e £ = u5t, and 5B = (dtB)St 
confirms equation ([Dip . 

As in WG can derive the MHD momentum equation 
using Hamilton's principle of stationary action. To cor- 
rectly incorporate the Lorentz force, we must take the 
Lagrangian to depend on B in addition to it, p, and s. 
Defining 



. = dC dC 
3 ~ du dB 

the general MHD Euler-Lagrange equation becomes 

dj 



(D4) 



dt 
-pV\ 



V • (uj) + Vu • j + 

(f) + i v » = < v * H > xB - (D5) 

Using the Lagrangian (|Lundgrenlll963|) . 

IBI 2 



C 



MHD 



u 



e(p,s) - 4> 



2m 



-, (D6) 



one may str aigh tforwardly check that the left-hand side 
of equation (|D5|) reproduces the Euler momentum equa- 
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tion, while the right-hand side of equation (|D5I) equals 
the Lorentz forcc0- 

Furthermore, for any Lagrangian depending on only u, 
B, p, and s, a conserved Hamiltonian, and generalized 
pressure tensor, follow respectively 



n = u-j - C, 
P 



dC dC \ D dC 



and satisfy 



d t H + V • (Hu + P.u) 



0. 



(D7) 
(D8) 

(D9) 



Thus, the Lagrangian derivation maintains energy con- 
servation. 

To eliminate the fast magneto-acoustic modes from the 
MHD equations, we first define a base state satisfying 
magnetohydrostatic balance, 

Vpo + *V* = (VX£ ° )XJ3 " . (D10) 
Mo 

In £14.11 we argued that sound waves rapidly equilibrate 
the pressure within a fluid parcel to that of its surround- 
ings. From this physical argument, we derived the GPI 
equations by constraining the system so p(p, s) = po(x). 
In the MHD context, magneto- acoustic waves rapidly 
equilibrate a fluid parcel's total pressure to that of its sur- 
roundings, where the total pressure comprises the sum of 
the gas pressure, and the magnetic pressure, 



U(p,s,B) = p(p, 



\B[ 

2p 



(Dll) 



Thus, we eliminate magneto- acoustic waves from the 
MHD equations by imposing the constraint 

C(p, s, B, x) = H{p, s, B) - n (x), (D12) 

and using the constrained Lagrangian. 

C = C M HB-XC(p,s,B,x), (D13) 

The full s et of equations associated follow from equa- 
tion ([Di~3|) . 



D P ^ 
Du — 

p— + {p- PQ )V4> = -v(rn A) + Avn 

B-V((l + A)fl) B VB 



Mo 



/'o 



DB 
~Dt 



BV • u = B • Vu, 



rn v-u + w vn = 

where 

] 

r = - 



(B ■ Vu) • B 



/'() 



JBJl 

2fj, p 



1 + 

2fj, p 



(D14) 

(D15) 
(D16) 
(D17) 

(D18) 



2 see un-published 2007 Mathematical Tripos, Part-III notes of 
Gordon Ogilvie for an excellent introduction to Hamiltonian MHD 



gives the average adiabatic exponent. If gas pressure 
dominates over magnetic pressure (2p ^> \B\ 2 /po), then 
r w Tj, whereas if magnetic pressure dominates over gas 
pressure (|B| 2 /^o ^ 2p), then T « 2. These equations 
conserve the same energy as the MHD equations, but 
have a modified pressure tensor, 



(n + rn A) i — (1 + \)bb. 



(D19) 



Impor tantly , we show in upcoming work that equa- 
tions (|D14[) - (|D18|) correctly reproduce the dynamics of 
Alfvcn waves, as well as two- and t hree-dimension al mag- 
netic buoyancy instability results (|Achesonll 19791) . 

Potential vorticity conservation no longer holds in the 
presence of magnetism. Recall from £13.31 that PV con- 
servation follows from the explicit invariance of the La- 
grangian with respect ot displacements in the form of 
equation (|80p . A major implication of PV conservation 
is that the two-dimensional pair {p, s} act as reduced 
phase-space coordinates for the three-dimensional phase- 
space momenta pu. Each reduction in phase-space de- 
pendence implies one conserved mome nta. 

In the case of magnetism, equation (|D3|) implies that 
the three-components of B also act as phase-space coor- 
dinates, bringing the apparent total to 5. It seems that 
magnetic field confuses the accounting of phase-space di- 
mension. However, 



V • B = 0, 



(D20) 



implies that only two of the additional components of 
equation (|D3|) vary independently. However, it still 
seems that we have one extra coordinate for three- 
dimensional dynamics. The final relation between the 
apparent four remaining coordinates comes from the con- 
servation of potential magnetism 



d_ 

dt 



VsB 



= 0. 



(D21) 



Equation (|D21[) contains no significant dynamical infor- 
mation, and simply constrains an over-specification of 
phase space accompanying equation (|D1[). No w however, 
equations ([S]) & © along with equation (|D21[) imply fully 
six-dimensional phase-space dynamics with no conserved 
momenta in the form of equation (|85|) . 

E. MODIFICATIONS TO THE MAESTRO 
CODE 

The numerical tests in ^required modifying the MAE- 
STRO code to implement the PI equations (|M f - (rilI)]) . 
The necessary modification follows from MAESTRO's 
implementation of th e velocity divergence con straint in 
the LM equations fsee lNonaka et al.ll2010ll2012l for more 
details) . 

The MAESTRO code assumes a velocity evolution of 
the form 



du 1 

ot p 



where the acceleration, 



(El) 



(E2) 
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Equation (|E1[) separates the accelerations due to pres- 
sure, and the accelerations due to inertia and gravity 
(and other sources in general). The MAESTRO code 
satisfies the constraint 



v-GSoti) = o 



(E3) 



using a divergence-cleaning scheme, evolving the velocity 
in two steps. 

To advance the velocity u n by At, MAESTRO first 
computes a™, from the solution at the nth time step, 
it'™', and /)". With this, the code advances the flow 
according to 



£(«+!) 



aAt. 



(E4) 

The code then removes (cle ans) the portion of u^ n+1 ^ 
not satisfying equation (|E3[) . This requires solving the 
elliptic equation, 



V- 



for x- The new flow solution 



u 



p 



(E5) 



(E6) 



satisfies equation (|E3|) . We can then identify x = p' /At. 

The PI equations also employ the constraint equa- 
tion (|E3|) . but use the velocity evolution equation 



du 
~dt 



Ay? (E. 

Jo 



(E7) 



The MAESTRO code easily implements this new velocity 
evolution equati on. F irst, we solve for u^ n+1 ' in the same 
way as equation (|E4[) . We then solve a modified equation 
for x, 



v.(fvx 

and update u^ n+1 ^ via 
..(«+!) _ 



;(«+!) 



fi (n+l) 



/?oVx 



(E8) 



(E9) 



We now identify x = p'/(/3oAt). To implement the PI 
equati ons inst ead of the LM equati ons, we repla ce equa- 
tions (jE5j) & dE6j) with equations flE8j| & (JEOJ), respec- 
tively. 

The distinction between the original MAESTRO im- 
plementation and our modifications amounts to simply 
replacing 



Vx -> A)V X . 



(E10) 



In both cases, we merely subtract a vector from the initial 
evolution equation (|E4[) . In both cases, the divergence 
of that v ector equals the divergence of the output from 
equation (|E4[) . The difference between the vortical com- 
ponents of these vectors provides the difference between 



the two methods. That is, even though 

v.(^vx) =v.(fv* 

it is also true that 



(Ell) 



unless 



V/3 x V X = 0, 



(E13) 



which is not satisfied in general. Therefore, we find that 
the PI and LM equations possesses non-trivial differences 
(see gBJ. 
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